% "Repricing Avalanche" Online Appendix D.1
% Obtain statistics of simulations of the two-productivity model
% with constant money growth policy
% July 2022
% Makoto Nirei

clearvars;
%load ../revision/da01/da01.mat;
load simu.mat;
dlogPm=zeros(maxmonth,pilength); L=zeros(maxmonth,pilength);    % monthly aggregate price increase
dlogPm(1,:)=sum(dlogP(1:calvom(1),:));                  % monthly price increase for the 1st month
idl=(nump(1:calvom(1),:)>0);                            % indicator for price rise
idlm=ones(calvom(1),pilength)-idl;                      % indicator for price cuts
cumcost_t = cumsum(nump.*(nump>1)*delta - nump.*(nump<-1)*delta_u); % cumulative menu costs
DL_t = cumcost_t -interp1([0; calvot],[zeros(1,pilength); cumcost_t],calvot-1/12,'previous'); % cumulated in the last month. 1st month will be discarded
init_menu=sum(calvot<1/12);
DL_t(1:init_menu,:)=ones(init_menu,1) * mean(DL_t(init_menu+1:end,:),1); % replace the first month by mean
mNpY=zeros(maxmonth-1,pilength);
mcons=mNpY; DL=mNpY; mCpY=mNpY; mlogm=mNpY; mY = mNpY; mw = mNpY; % pre-allocation
pc1=(1-alpha_U*(1-gamma_U))/(1-alpha_U)/(1-gamma_U);
m_ss=zeros(1,pilength); c_ss=m_ss; n_ss=m_ss;
for i_pi=1:pilength
    ss=ss_pi{i_pi};
    m_ss(1,i_pi)=ss.m;
    c_ss(1,i_pi)=ss.C;
    n_ss(1,i_pi)=ss.N;
    menu_ss(1,i_pi)=ss.MenuCost;
end
mlogm(1,:)=log(m_ss(1,:));

options = optimset('Display','off');
for i=2: maxmonth
    dlogPm(i,:)=sum(dlogP(calvom(i-1)+1:calvom(i), :)); % monthly price increase
    mlogm(i,:)= mlogm(i-1,:)-dlogPm(i,:)+ piv*month_time; % monthly real money balances
    idl=(nump(calvom(i-1)+1:calvom(i),:)>0);            % indicator for price rise
    idlm=ones(calvom(i)-calvom(i-1),pilength)-idl;      % indicator for price cut
    DL(i-1,:)=DL_t(calvom(i-1)+1,:);
    mNpY(i-1,:) = (sum(NpY(calvom(i-1)+1:calvom(i),:).*diff(calvot(calvom(i-1)+1:calvom(i)+1)))...
        +NpY(calvom(i-1),:)*(calvot(calvom(i-1)+1)-month_time*(i-1))...
        -NpY(calvom(i),:)*(calvot(calvom(i)+1)-month_time*i))/month_time; % monthly averaged NpY
    mCpY(i-1,:)= 1-DL(i-1,:)/n/month_time;

    compu_Y = @(x) rho*alpha_U.* exp(mlogm(i-1,:)) .* ( (mCpY(i-1,:) .* x).^alpha_U...
        .* (1- mNpY(i-1,:).* x).^(1-alpha_U) ).^(1-gamma_U) - mCpY(i-1,:).*x;
    mY(i-1,:)=fsolve(compu_Y, c_ss, options);
    mw(i-1,:)= mCpY(i-1,:).*mY(i-1,:) *(1-alpha_U)/alpha_U ./ (1-mNpY(i-1,:).*mY(i-1,:));
end

i_pi=2;
num_2years= floor((maxmonth-init_menu)/24);
std_mw2=std(log(reshape(mw(init_menu:init_menu-1+24*num_2years,i_pi), 24, num_2years)));
[mean(std_mw2) std(std_mw2)]

num_5years= floor((maxmonth-init_menu)/60);
std_mw5=std(log(reshape(mw(init_menu:init_menu-1+60*num_5years, i_pi), 60, num_5years)));
[mean(std_mw5) std(std_mw5)]

% dlogPy=zeros(maxmonth-11,pilength);
% for i=1:maxmonth-11
%     dlogPy(i,:)=sum(dlogPm(i:i+11,:));                  % monthly, year-on-year price increase
% end
% stdP = std(dlogPy);                                     % volatility of annual inflation
% meanP = mean(dlogPy);                                   % average of annual inflation
% 
% % stats for various pi
% disp([piv; meanP; stdP; negnump]);
% 
% disp('std of log(c) for various pi');
% disp(std(log(mcons)));
% 
% disp('std of monthly log(N/Y)');
% disp(std(log(mNpY)));
% 
% disp('mean, std, min, max of monthly d(log(P))');
% disp([mean(dlogPm(:,i_pi)) std(dlogPm(:,i_pi)) min(dlogPm(:,i_pi)) max(dlogPm(:,i_pi))]);
% 
%i_pi=2;

len=length(dlogPm(:,i_pi));
figure(20),
set(gca,'FontSize',20);
yyaxis left
%plot(100:len,cumsum(dlogPm(100:len,i_pi))./(((0:len-100)'*piv(i_pi)*month_time) + sum(dlogPm(1:100,i_pi))),'-');
plot(1:len,cumsum(dlogPm(1:len,i_pi)) - (0:len-1)'*piv(i_pi)*month_time,'-');
yyaxis right
plot(2:len,mw(:,i_pi),'--');
xlabel('Months');
legend('log(P)-pi*t','real wage');
title('Typical paths of monthly log-price and real wage');

% figure(21),
% set(gca,'FontSize',20);
% stairs(calvot(100:200)-calvot(100),cumsum(dlogP(100:200,i_pi))+1);
% xlabel('Time (unit time is a year)');
% ylabel('log P');
% xlim tight;
% title('Typical path of log P(t)');
% 
% figure(22),
% set(gca,'FontSize',20);
% plot(0:100,dlogPm(100:200,i_pi));
% xlabel('Months');
% ylabel('logP(t+1)-logP(t)');
% title('Typical path of monthly inflation');
% 
% corrcoef(dlogP(100:end-1,2),dlogP(101:end,2))
% corrcoef(dlogPm(100:end-1,2),dlogPm(101:end,2))